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We study solutions of the one-dimensional Gross-Pitaevskii equation to better understand dynam¬ 
ical instabilities occurring in flowing atomic condensates. Whereas transonic stationary flows can 
be fully described in simple terms, time-dependent flows exhibit a wide variety of behaviors. When 
the sound speed is crossed once, we observe that flows analogous to black holes obey something 
similar to the so-called no hair theorem since their late time profile is stationary and uniquely fixed 
by parameters entering the Hamiltonian and conserved quantities. For flows analogous to white 
holes, at late time one finds a macroscopic undulation in the supersonic side which has either a fixed 
amplitude, or a widely varying one signaling a quasi periodic emission of solitons on the subsonic 
side. When considering flows which cross the sound speed twice, we observe various scenarios which 
can be understood from the above behaviors, and from the hierarchy of the growth rates of the 
dynamical instabilities characterizing such flows. 
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I. INTRODUCTION 

A recent experiment [1] performed with a flowing ultra cold atomic Bose condensate has observed several features 
which are in agreement with predictions concerning the so-called black hole laser instability. This instability was 
found in Ref. [2] and further studied in Refs. [3-6]. It occurs when the speed of a one-dimensional steadily flowing 
condensate is supersonic in a finite internal region surrounded by two subsonic external regions. The internal region 
acts as a resonant cavity as it contains a discrete set of trapped negative-energy phonon modes. In addition, the 
mixing of these cavity modes with the standard positive energy modes which propagate in the external regions leads 
to a self-amplified super-radiance. As a result, the frequency of the trapped modes acquires a nonvanishing imaginary 
part that fixes the growth rate of their amplitudes. Finally, it should be pointed out that the mode mixing occurring 
near the two sonic horizons limiting the supersonic region is the analog version of the Hawking effect [7, 8]. 

Several important observations have been made in Ref. [1]. In the internal region, the late time evolution of g 2 , 
the density-density correlation function, is clearly governed by a single complex-frequency mode. This is in accord 
with Ref. [4] where it was predicted that the mode with the highest growth rate should dominate the late time 
behavior of the system. In addition, the spatial properties of 52 exhibit both the fixed nodes of the trapped mode, 
and correlations with the emitted phonons of positive energy. These two observations are in good agreement with the 
theoretical analysis of Ref. [5] . Interestingly, Steinhauer also observed that the ensemble averaged density develops a 
growing pattern of fluctuations that has essentially the same spatial profile as that of g 2 - This observation is a surprise 
because the linear treatment of Ref. [5] predicts that the mean value of the density fluctuations should vanish at all 
times, unless these are seeded by classical perturbations present in the initial conditions. In fact, at linear order, the 
growth of the mean value from zero is forbidden by a Z 2 symmetry relating solutions of the Bogoliubov-de Gennes 
equation. An alternative explanation of this observation, which does not refer to inhomogeneous initial conditions, 
arises from the nonlinear effects associated with the exponential growth of fluctuations. As conjectured in Ref. [9], 
and verified in Ref. [6], the Z 2 symmetry breaks down when nonlinear effects become significant. As a result, even 
when it initially vanished, the mean value should grow exponentially in time in the situation of Ref. [Ij. 

In this paper, we aim to clarify the situation and complete former studies. The novelty consists in analyzing the 
temporal evolution of density perturbations in transonic flowing condensates. The parameters describing our systems, 
such as the inter-atomic coupling, are all constant in time. As a result, the time dependence of the solutions is seeded 
by initial conditions. Yet, our analysis reveals the nontrivial aspects of the nonlinear dynamics of transonic flows. To 
identify the key elements, we first analyze monotonic transonic flows which cross the sound speed only once. There are 
two such cases, corresponding to flows which accelerate or decelerate in the direction of their velocity. When the flow is 
accelerating (decelerating), one gets a black (white) hole flow in that the wave number of incoming counterpropagating 
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waves is reduced (increased) when crossing the sonic horizon [7]. Even though the set of stationary states is the same 
for black and white flows, we find that their time evolutions are very different. 

When considering black hole flows, at early time we observe an emission of dispersive shock waves [10]. The resulting 
flow is stationary and asymptotically homogeneous on both sides. These two properties, along with conservation of 
the total mass and energy, uniquely determine the late-time flow. This configuration thus seems to act as an attractor 
in the space of solutions. This is very reminiscent of the gravitational black hole “no hair theorem” [11], which states 
that stationary black holes in (3-1-1) dimensions are uniquely characterized by their mass, angular momentum, and 
electric charge. For white hole flows we found two different behaviors which are related to the above mentioned Z 2 
symmetry. Depending of the sign of the detuning with respect to the homogeneous solution, white hole horizons 
emit an undulation in the supersonic region which either has a fixed amplitude, or which widely varies because it is 
accompanied by a seemingly infinite number of solitons in the subsonic region. 

We then study flows which twice cross the sound speed and which are subsonic in the external regions. By a 
linear stability analysis, we determine the stability level of the relevant set of stationary solutions. We show that the 
growth rate of the most unstable mode introduces a clear hierarchy in this set. Combined with the above analysis of 
monotonic flows, this hierarchy allows us to understand the generic properties of time-dependent solutions. Particular 
attention is accorded to the breakdown of the Z 2 symmetry as it indicates when nonlinear effects can no longer be 
neglected. 

To complete the analysis undertaken in Refs. [5, 6], we also study the consequences of a small detuning, i.e., a small 
change of the system parameters with respect to those characterizing stationary homogeneous flows. Even though the 
linear series of stationary solutions is affected by a (small) detuning, their physical properties and their associated set 
of dynamically unstable modes, are both mildly affected. In other words, the properties of solutions obtained with 
fine-tuned parameters are generic in character. 

The appendices give complementary information on various aspects of this work. Appendix A shows the conse¬ 
quences of breaking the Z 2 symmetry on the density-density correlation function in flows with two horizons. Appendix 
B details the calculations giving the stationary solutions and dynamically unstable modes in such flows. Finally, Ap¬ 
pendix C gives the main properties of the dispersive shock waves emitted by black hole flows. 


II. TIME EVOLUTION OF A SINGLE BLACK OR WHITE HOLE FLOW 

This section is organized as follows. We first briefly review the main properties of the stationary solutions which 
are transonic and asymptotically bounded. Interestingly, the complete set of such solutions is characterized by two 
parameters related to the amplitude of the density modulations in each asymptotic region. This naturally leads to the 
notion of fine-tuned solution, for which both of these quantities vanish. Using numerical simulations, we then explore 
the time evolution when starting with initial conditions which do not coincide with a stationary solution. We show 
that white hole and black hole flows, although symmetric in time independent cases, behave very differently in time. 


A. Tuned and detuned stationary solutions 

We consider a flowing atomic Bose-Einstein condensate in I-l-I dimensions. We work in the so-called quasicondensate 
regime [12] where it is legitimate to locally treat the wave function as a classical held, although in a strict sense no 
Bose-Einstein condensation occurs. In this regime, the condensate wave function ip satishes the Gross-Pitaevskii 
equation (GPE) 

idtip =+ V(x)i> + (I) 

We work in a system of units in which the reduced Planck constant h and the atomic mass are equal to 1. In this 
system, the local value of the healing length is ^{x) = {\'ip{x)\y/2g{x))~^. The unit length is then dehned by imposing 
that \ip\ —>■ I for X —>■ — 00 . We look for stationary solutions of the form 

V-(x,t) = e-“‘/(x)e*®("), (2) 

where w S K and /, 9 are two real-valued functions. The GPE then gives 

dlf = 2gf-2gf + ^, 
where g,{x) = uj — V(x) and where J = f^dxO is the conserved current. 


( 3 ) 
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FIG. 1. (Color online) Solutions of Eq. (4) for J = -^8/3, fb = \/2, and fp = 1. The two solid lines give the homogeneous 
supersonic solution f = fp (red, middle), and the subsonic one f = fb (blue, top), obtained respectively for C = Gmin = 19/3 
and C = Gmax = 20/3. The other solutions are obtained for C = 6.4 (dashed, orange), C = 6.6 (dotted, green), and 
C = 6.666666 (dot-dashed, cyan). One sees that the solutions near Gmin can be seen as adding a periodic undulation on top 
of the supersonic flow f = fp, whereas those near Gmax can be seen as adding a train of deep solitons on top the subsonic flow 
/ = /{,. One also sees that there is a smooth transition from super to subsonic flows. 


When fj,,g are constant, and when there exist two homogeneous positive solutions d^f = 0 of Eq. (3). 

One of them, which we shall call fp, corresponds to a supersonic flow, while the second one, fb, corresponds to a 
subsonic flow [10]. Integrating Eq. (3) gives 

= + (4) 

where G S K is an integration constant. Bounded solutions exist provided Gmin < G < Gmax, where the two extremal 
values are given by setting dxf = 0 in Eq. (4): Gmin corresponds to / = fp, and Gmax to / = fb- For small 
G — Gmin > 0, the solution contains a small-amplitude sinusoidal undulation on top the homogeneous solution / = fp, 

with a wave vector ko = 2y^J2//4 — When increasing G, the wavelength increases and the solution is deformed; 

see Fig. 1. For G —>■ Gmax from below, one obtains widely spaced solitons, which have a very low density at their core 
and which are separated by a region in which f ^ fb- 

By increasing G from Gmin to Gmax, there is a smooth transition from the homogeneous supersonic solution to the 
subsonic one without any clear separation. A possible criterion to define a supersonic regime, and a subsonic one, is 
to require that |u| /c > 1 (supersonic) or ju] /c < 1 (subsonic) over more than a fraction 2/3 of the domain, where 
V = J/p \s the local fluid velocity, and c = ^gf is the local sound speed. This choice defines two new critical values 
of G which respectively give the highest value of G for supersonic flows Gm£f'' and the lowest one for the subsonic 
regime 

In what follows we restrict our attention to transonic stationary flows which are engendered by piecewise constant 
/i and g > 0 having both a single discontinuity at a: = 0. We shall denote with an index quantities evaluated 
for a; > 0 and with an index ” quantities evaluated for a; < 0. We shall also assume that the flow is supersonic on 
the right side. Then, for J > 0 (J < 0), we get a black (white) hole flow. For more details about the correspondence 
between transonic flows in BEG and black holes we refer to Refs [7, 8]. 

In this article, we call a set of 4 -|- 2 parameters V±,g±, J, and w “fine-tuned” if the densities fb - and fp^+ obey 

fb,- = fp,+ , (5) 

which means that there exists a uniform density solution which is subsonic for a; < 0 and supersonic for a; > 0. 
These peculiar sets of parameters have been used in many works to study the analog Hawking radiation and related 
effects; see [6, 8, 13-17]. It should be noticed that the system parameters V±,g± are truly constant and appear in the 
Hamiltonian, whereas the values of J and ui will in general vary in time when considering nonstationary solutions. 

In fact, for each set of the 4 parameters V±,g± such that V+ — V- and g+ — g- have opposite signs, there exists a 
linear series of solutions where the current J can take any value and where the frequency w of the solution is found 
by requiring that Eq. (3) be satisfied in one of the two regions. In the steep regime limit we are using, it is then 









4 



FIG. 2. (Color online) The spatial profile of the square root / = |^/)| of the mean density as a function of x for two stationary 
solutions with the fine-tuned parameters J = -y/S/S, g+ = 1, g- = 8, = 7/3, and fi- = 28/3 (solid, black), and when 

imposing that the solution contains no soliton in the subsonic region. The amplitude of the undulation in the supersonic side 
is fixed by the integration constant of Eq. (4) for a; > 0, C+ = 6.4. The dashed blue curve represents the shadow soliton (left) 
or the soliton (right) which coincides with the solution for a; < 0. Notice that the density of the first one is larger than that of 
the homogeneous solution, thereby reducing the local speed of the flow. 


automatically satisfied in the other one. In addition, in this regime, irrespectively of J, / is given by 


fi.t. 


IV- - V. 


+ 


ff+ - 9- 


( 6 ) 


When abandoning the steep regime limit, and considering V and g given by smooth and monotonic profiles, Eqs. 5 
and 6 will no longer apply. But we strongly conjecture that there will still be a one-parameter family (labeled by 
J) of stationary solutions which are asymptotically uniform on both sides. An interesting example is the water fall 
configuration discussed in Ref [18]. We are currently trying to prove this conjecture. ^ In this perspective, Eqs. 5 
and 6 should be conceived as particular expressions of this general fact in the steep regime limit. 

Given some fine-tuned parameters in this limit, it is interesting to study the set of stationary transonic flows which 
are asymptotically bounded on both sides. This set is of dimension two, as it was the case when working with uniform 
V and g. In that case, stationary solutions were characterized by the constant C up to translations in x. When V 
and g vary, the set is more appropriately parameterized by the two integration constants C± of Eq. (4) evaluated on 
each side of the discontinuity of V and g. (Notice that for each choice of (<7+, C-), the number of stationary solutions 
is generally larger than 1, as the corresponding lines in the phase portrait in general cross each other several times; 
see Fig. 16. In this set, we will mostly focus on solutions which do not reach the bottom of the soliton in the subsonic 
region.) 

When further imposing that the solution is asymptotically constant on the subsonic (—) side, one must choose 
C- = Then C'_|_ fixes the amplitude of the periodic density modulations on the supersonic side. At the 

linear level, in the supersonic region the undulation is a sinusoid with wave vector fcg = — c\. In the subsonic 

region, it is exponentially decaying as with kq = The flow velocities and the sound speeds in each 

region are respectively given by = J//p +, V- = J/f^_, c+ = y/g+f^, and c_ = y/gZfj^. The condition that 
the perturbation is bounded in the subsonic region fixes, up to a sign, the density perturbation 6f = f{x) — 1 in the 
supersonic region; see Fig. 2. 

At the linear level, there is thus a Z 2 symmetry between these two solutions with the same values of C±. At the 
non linear level, this Z 2 symmetry is broken in the following sense. The solutions for which f{x = 0) > fh,- have 
a fraction of a shadow soliton for a; < 0, while solutions for which f{x = 0) < fb- have a fraction of a soliton. 
Importantly, the former has a lower energy than the later; see Ref. [6]. [A proper evaluation of the energy can be 


^ We numerically verified it for couples V,g related by the condition V{x) + g(x)fQ = F, where fg and F are constant. We also studied 
cases where V and g are given by sums of constants and hyperbolic tangents with slightly different slopes so that the above condition 
is not exactly satisfied. For each J, we always found one solution which is asymptotically homogeneous on both sides, with asymptotic 
densities differing by a term of the order of the relative difference between the two slopes. 

^ In addition to this choice, there exists a continuous class of solutions which contain soliton trains for Cmin,— < C— < _, where 

the upper bound is the value above which the flow for a; < 0 is no longer considered subsonic. 
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done when including a second horizon, so as to have a finite supersonic region. Notice also that here is a third series 
of solutions which have a larger fraction (more than 50%) of a soliton in the region x < 0. However, as these solutions 
have a larger energy, and are not directly connected to linear perturbations around the homogeneous configuration, 
we shall discard them.] 

Let us now briefly present the (small) modifications of the set of solutions when working with detuned parameters, 
for which Eq. (5) does not hold. To start the analysis, we assume that the detuning is small. This means that the 
above six parameters are such that a nearly uniform transonic solution exists. 


lA- - fp,+l < fp-,fp,+, \A- - A +\> I/p - - /p.+l • 


( 7 ) 


We used the values of / rather than the 6 parameters because the condition to be fine-tuned involves only and 
fp^+; see Eq. (5). In the following, it will be convenient to refer to the “sign” of the detuning in the following sense. 
We call the detuning positive if — /p__|_ > 0 and negative otherwise. 

The novelty introduced by the detuning is that the amplitude of the undulation in the supersonic region cannot be 
set to zero [16]. Its nonvanishing minimum value can be correctly conceived as forced by the detuning. To first order 
in \fb,- — fp,+ \, the minimum amplitude of the undulation in the supersonic region is equal to 


^min 


\A- 


/p.+l 



+ o 


f (A--fp,+f ] 
Y fb,- + fp,+ j 


( 8 ) 


As the next sections and subsections heavily rely on nonlinear solutions of the GPE, it is of interest to see how 
to characterize the minimum-amplitude undulation at the nonlinear level. We find the solutions are qualitatively 
unchanged provided the two following conditions are satisfied: 


0 < 


M-e - 1 ^- 
9+-9- 


< fh 


and 


9^ fA , - 7 ^ , iP+ - 9-Y ^ 9+ 
2 


A 


< w/6,+ + 


fb,- ' 2{g+-g-) 


fl+ 


( 9 ) 


( 10 ) 


A straightforward calculation shows that the value of (7+ giving an undulation with the smallest possible amplitude 
is 

\2 

(11) 


r = A+-p-y 

2 /b.- ^ 2 _ 2ig+-g-)- 


The important result to retain from this subsection is that the space of stationary transonic solutions is of dimension 
two. The two parameters characterize the amplitude of the density perturbations on each side. In addition, when 
fb,- fp,+ , the asymptotically uniform solution for x —> —oo necessarily has a nonvanishing undulation for a; > 0. 

Notice that stationary solutions in / depend on J only through A. Therefore, in this subsection the sign of J plays 
no role. However, when considering the time evolution of a transonic flow, this sign affects the dynamics in a crucial 
way. Indeed, changing the sign of J changes the sign of the group velocity of the undulation; to wit, its velocity is 
oriented towards the horizon for a black hole flow (J > 0) and away from it for white holes (J < 0). As a result, 
white hole flows can and will emit undulations, whereas black hole flows cannot. 


B. Black hole flow: An analog no hair theorem? 

Since black hole horizons cannot generate a stationary solution by emitting an undulation, it is not straightforward 
to guess what will be the end point of the evolution when starting with some arbitrary initial condition. To determine 
this time evolution, we numerically solved the GPE for detuned sets of parameters fp,+ fb,-- 

The GPE was integrated on a torus of length SOtt with periodic boundary conditions, with a space step of 0.003 
and a time step of 0.002. In order to check the stability of our code, we first chose uniform values of p, and g and took 
initial conditions corresponding to solitons with known velocities. We checked that their propagation was consistent 
with analytical predictions and saw no deformation for times larger than the duration of the simulations we report 
below. We also checked that the code gives the correct stationary solutions when using the settings of Ref. [6] . 

To simplify the analysis, the initial condition of the simulations detailed in the forthcoming plots is taken to be 
homogeneous and obeying f(t = 0) = fb,-- We also ran simulations with f{t = 0) = fp,+ and found very similar 
results. In Fig. 3, we clearly see that three shock waves are emitted. The main properties of the emitted dispersive 
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FIG. 3. (Color online) The spatial profile / = |i/)| as a function of x evaluated at f = 0 (blue, solid), t = 2 (orange, dashed), 
t = 10 (green, dotted), and t = 30 (red, dot-dashed) for a detuned black hole flow with J = •y/8/3, ft,- = 1, fp,- = 0.7, 
= 1.5, and /p,+ = 1.2. The discontinuity of V and g is located at a: = 0. We start at t = 0 from a homogeneous 
configuration f — 1 and d^O — J > 0. We observe that the final value of / at a; = 0 (on the sonic horizon) is / ~ 1.03 in very 

good agreement with the fine-tuned value /f.t. = \/ 2930900 ~ 1-0305 given by Eq. (6). 


shock waves are given in Appendix. C. It is shown that, to linear order, shock waves propagate with velocities t) ± c 
and are followed at late times in the supersonic region by an oscillation tail with wave vector kc where the dispersion 
relation has a horizontal tangent; see Eq. (24) in Ref. [7]. This can be understood as follows. To linear order, the initial 
configuration may be seen as a perturbation with respect to the exact homogeneous solution. This perturbation is a 
superposition of modes with real frequencies and wave vectors. When sending t —>■ 00 , only the modes with a vanishing 
group velocity will not have propagated away from the horizon. Close to the horizon at a; = 0, their contribution to 
the density perturbation decays as see Appendix C. This result is in accord with our simulations, where we 

verified that the properties of the oscillations in the tail of the second shock wave emitted to the right are independent 
of the detuning parameter /{,__ — fp,+ . 

Interestingly, for all simulations, we observed that the emission of dispersive shock waves is always followed by a 
stationary profile f{x) which is homogeneous across the horizon. A straightforward calculation using Eq. (3) and the 
definition of fi shows that the final value of / is given by /t.t. of Eq. (6). In other words, the temporal evolution 
of the GPE sends the flow towards the unique stationary solution which is asymptotically homogeneous on both 
sides, and has a current compatible with the conserved quantities of the GPE; see Appendix C.^ It thus seems that 
one-dimensional analog black holes lose their hair. 

We now give a qualitative argument which validates these numerical observations. We assume that the solution for 
t f 00 

1. remains transonic, with v < c for a; < 0 and v > c for a: > 0, 

2. becomes stationary dtf ^ 0, 

3. preserves the sign of the group velocity of the undulation, and 

4. preserves the inequality fb^+ < fp,+ - 

From the hrst two assumptions, the solution in the region a; > 0 is a (nonlinear) superposition of a homogeneous one 
and an undulation. From the third one, the group velocity of the undulation is negative, so it can not be produced 
at a; = 0, and its amplitude is necessarily 0.^ So, / must be exactly uniform in the region a: > 0 at late times. When 


® Using the parameters of Fig. 3, we checked that this property remains true with smooth functions V —^x^±oo V± and g —>^->±00 g±- 
As in the steep horizon limit, we found that the late-time configuration (reached again after the emission of three shock waves) always 
coincides with one of the stationary solutions with asymptotically uniform densities at a; —> ± 00 . 

^ This is true unless the undulation is produced by a shock wave previously emitted. At linear order, this is not allowed. Thus this 
possibility is a loophole of the nonlinear version of our present argument. 











7 


f 



FIG. 4. (Color online) The density profile / = \tp\ as a function of a; at t = 5 (blue, solid) and t = 30 (orange, dashed) for a 
detuned white hole configurations with J{t = 0) = —^/8/3, fb,- = 1, fp,- ~ 0.7, fb,+ = 1.5, and /p,+ = 0.9, when the detuning 
sends the solution towards the shadow soliton. We start at f = 0 from a homogeneous configuration / = 1 and 0{x) = Jx. At 
early time, we clearly see the small shock wave emitted in the subsonic left region, and the growth of the undulation. At late 
time, we observe that the undulation amplitude saturates and that its nodes are fixed. 


applying the matching conditions at a; = 0, we find two possibilities: / is either uniform also in the region x < 0 
or contains half a soliton connecting the region x —oo to x > 0. The second possibility is forbidden by the last 
assumption. So, / must be uniform for t —>■ oo. In a future work, we hope to be able to give a more rigorous proof 
which would provide an analogous black hole no hair theorem, valid at the nonlinear level as based on solutions of 
the GPE. 


C. White hole flows 

We now consider transonic flows with J < 0. A linear analysis reveals that the undulation propagates away from 
the sonic horizon in the supersonic region and that its amplitude has the tendency to grow in time [9, 14]. For 
fine-tuned white hole flows, it was shown that such an undulation is produced by an infrared instability and that the 
saturation mechanism involves the suppression of the instability by the growth of the undulation [16]. For generic 
initial conditions, we thus expect that an undulation with a large amplitude will be emitted. 

Figure. 4 shows results for simulations with a positive detuning parameter — /p^+ = 0.1 such that the initial 
condition sends the solution towards the shadow soliton. We set J = --^8/3 at t = 0, fb - = 1, fp - = 0.7, /(,_+ = 1.5, 
and /p_+ = 0.9. At < = 0, the solution has a homogeneous density with / = 1. At early times it shows a damped 
undulation emitted from the white hole horizon in the supersonic region x > 0, as well as a small dispersive shock 
wave emitted in the subsonic region x < 0. As expected from the linear analysis, the amplitude of the undulation 
grows in time. For large values of x, its amplitude seems to be proportional to if jx^. We observe that it saturates 
with a relative amplitude of « 28% of fb,-- It should also be noticed that while the right front of the undulation 
propagates to the right, its nodes are fixed. Close to x = 0, this undulation has a crest, so the solution for x < 0 is 
part of a shadow soliton; see Fig. 2. 

Figure. 5 shows results for simulations with an opposite detuning, namely fb,- — fp,+ = —0.2. The parameters 
are the same as above, except that /p^_|_ = 1.2. At early times, the time evolution is similar to that of the previous 
case. However, the late-time regime is very different, as a seemingly infinite series of (equally spaced) solitons 
which propagate in the subsonic region is emitted from the white hole horizon. Other periodic trains of solitons are 
superposed on this series, with a much lower amplitude. We also observed that the amplitude of the undulation to 
the right is significantly affected by the successive emissions. Contrary to the solution of Fig. 4, the present solution 
is apparently not stationary at late time. Similar soliton trains have been thoroughly studied in Refs [19, 20]. 

We leave a precise study of these solutions to a future work. The important point to retain from this subsection is 
that, although the short-time dynamics of the two cases are similar, saturation effects introduce qualitative differences. 
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FIG. 5. (Color online) The density profile / = \ip\ as a function of a; at t = 4 (blue, solid) and t = 30 (orange, dashed) for a 
detuned white hole configuration with the same initial condition and the same parameters as in Fig. 4, save for the value of 
fp,+, which is now equal to 1.2. At late time, on the left side, one sees superposed trains of equally spaced solitons emitted 
periodically. The largest ones are very deep, as the density in their core is about 10% of the mean density. (The three solitons 
have the same depth. The apparent differences are due to the poor resolution of the plot.) The small ones are not clearly visible 
as their amplitude is small, of order 10%. On the right side, the amplitude of the undulation widely varies, as a consequence 
of these emissions. 


In one case, they simply stop the growth of the undulation, falling on a locally stationary solution. The corresponding 
solution in the region a; < 0 is a fraction of a shadow soliton. In the other case, we observed a seemingly infinite 
number of solitons. This discrepancy is a dynamical consequence of the breaking of the Z 2 symmetry which is present 
when dealing with solutions of the Bogoliubov-de Gennes equation, as discussed in the former section. When studying 
the dynamics of a black hole laser, we shall observe other manifestations of this symmetry breaking. 


III. BLACK HOLE LASER INSTABILITY 

We study black hole laser configurations in which fi and g are piecewise constant with two discontinuities, and take 
the same values at x —>■ 00 and x —>■ — 00 . We aim to complete the analysis undertaken in Ref. [6]. In that work it 
was shown that the homogeneous (fine-tuned) stationary solution possesses a discrete set of unstable modes which 
grows when increasing the interhorizon distance. In addition to the solution with homogeneous density which exists 
for any positive value of L, four inhomogeneous “connected solutions” were studied. These can be made arbitrarily 
close to the homogeneous one by varying L, and are of 4 types, called “sh-sh’,’ “sol-sol,” “sh-sol,” and “sol-sh.” Their 
names reflect their behaviors for x < —L and x > L, i.e., whether they contain part of a shadow soliton (“sh”) or of 
a soliton (“sol”) at each horizon; see Fig. 2. Here, we determine their instability level, and establish that there is a 
clear hierarchy governed by the growth rate of the most unstable mode. 


A. Linear dynamical stability 

We work with fine-tuned parameters and consider a setup with two discontinuities in V and g, located at x = ±L. 
For simplicity, we assume that V(x < —L) = V{x> L) and g{x < —L) = g{x > L). We denote quantities evaluated 
in the internal region —L < x < L by an index “int” and quantities evaluated in the regions |x| > L by an index 
“ext.” To study the dynamical stability of a stationary solution, we look for the set of asymptotically bounded modes 
(ABM), solutions of the linearized wave equation, whose angular frequencies A = w -b iF, (a;,F) S have positive 
imaginary parts F > 0. These modes thus grow exponentially in time, triggering a laser effect [4]. 

We note that in Ref. [6] two different types of unstable sectors were found. The usual one, called “nondegenerate”, 
with Lo and F both nonvanishing, corresponds to a complex unstable harmonic oscillator. We observed that this case 
was preceded (when increasing L) by a real oscillator with a; = 0, called “degenerate.” When the background flow 
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FIG. 6. Energy E and dynamical stability of the first nonlinear solutions found when increasing the interhorizon distance 2L. 
We work with “fine-tuned” black hole laser configurations with J = i/8/3, fb,mt = %/2, /p,ext = l/%/2, and /(,,ext = /p,int = 1- 
The degree of instability is indicated by the style of the curve: dotted indicates stable, small dashes indicate one degenerate 
dynamical instability, medium dashes indicate one nondegenerate dynamical instability, large dashes indicate one degenerate 
and one nondegenerate instabilities, and continuous indicates two nondegenerate dynamical instabilities. The insert shows a 
zoom on the point where the second “sol-sol” solution appears. As could have been expected, when a new instability occurs 
for increasing L, the degree of instability of the homogeneous solution is transmitted to a new inhomogeneous solution with a 
smaller energy. Therefore the ground state is the “sh-sh” solution with n = 1. 


is described by the fine-tuned homogeneous solution, n unstable modes were found for L 2 n -2 < L < L 2 n, where 
the values of for n S N are recalled in Eq. (B13), and by convention Ln — 0 for n < 0. Moreover, one of these 
modes corresponds to a degenerate instability if ^ 2^-2 < L < L 2 n-i, while they all correspond to nondegenerate ones 
otherwise. In the present work, we extend the analysis to inhomogeneous solutions which are smoothly connected to 
the homogeneous one. In the body of the text we only present the main results. The details of the analysis can be 
found in Appendix B 2. 

Figure 6 shows the energy E (defined in Appendix B 1) of the first connected stationary solutions which are 
homogeneous in both asymptotic regions, along with their degree of instability. We notice that, for each n S N, there 
exists a series of solutions for L > which has the same set of ABM (same numbers of degenerate and nondegenerate 
ones) than the homogeneous solution for L < Ln. This solution merges with the homogeneous one for L = L„. So, 
each time the degree of instability of the homogeneous solution is increased, a new inhomogeneous solution which 
preserves the number and type of dynamical instabilities continuously emerges at L = As a result, for any value 
of L there is only one dynamically stable inhomogeneous solution. It corresponds to the “sh-sh” solution with n = 1. 
This result remains valid when including the stationary solutions which are not connected to the homogeneous one, 
as all these solutions are dynamically unstable. 

To complete the analysis, we now study the relative magnitude of the instability of the above flows. Figure 7 shows 
the imaginary part of the frequency of the most unstable modes on the homogeneous solution, as well as on “sh-sh,” 
“sh-sol,” and “sol-sh” solutions. It can be seen that for a fixed value of n, the instabilities on “sh-sh” solutions are 
much milder than those on the homogeneous solution, while those on “sh-sol” and “sol-sh” solutions are stronger. 
Consequently, for sufficiently short time scales the “sh-sh” solutions can be seen as stable, while the other solutions 
are strongly unstable. The interested reader can find more details in Appendix B 2. The implications of this hierarchy 
will become clear when studying time-dependent effects. 
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FIG. 7. As a function of half the inter horizon distance L, we represent the imaginary part of the frequency of the most unstable 
mode on the homogeneous solution (circles), the first sh-sol solution n = 1 (empty squares), the second sh-sol solution with 
n = 2 (filled squares), the sh-sh solution with n = 2 (empty triangles), and the sh-sh solution with n = 3 (filled triangles). One 
clearly sees that sh-sh solutions are only mildly unstable compared to the other ones. 


B. Stationary solutions in the detuned case 


Before studying time-dependent effects, it is worth verifying that a small detuning does not significantly affect the 
main conclusions of the above analysis. To start, we remind the reader that the stationary solutions of Eq. (3) can 
be written in terms of elliptic Weierstrass functions. We refer to Ref. [6] for details on how these solutions can be 
obtained. The same method works, with minor modifications, in the detuned case fb,ext ^ fp.'mt, see Appendix B 1. 

The energy Eq. (B8) of the first solutions is shown in Eig. 8. Eor definiteness, we assume the detuning is small, in 
the sense that 




p,ini\ f 2 


<r fh <r 


( 12 ) 


In other words, the subsonic homogeneous density in the external region is between the two extremal densities for the 
stationary soliton in the internal region. We also impose that 


Jp,ext 


f? 


p,ext 


+ /f 


< fp, 


p,in • 


6,ext 


(13) 


These two conditions are always satisfied provided the parameters are sufficiently close to being fine-tuned. Then, the 
main difference with respect to the fine-tuned case is the following: When fb^ext = /p.int, the homogeneous solution 
/ = /h.ext exists for all values of L and is connected to an infinite number of series of solutions. When fb,ext 7^ /p,int, 
and when increasing L, each series of solutions is now connected only to a finite number of other series, as can be 
seen in Eig. 8 for the first few solutions. We checked this remains true when including all solutions, as can be easily 
deduced from plots of the phase portrait of Eq. (3); see Eig. 16. Continuity of the set of solutions in the limit of 
fine-tuning is recovered when noticing that, for very small detunings, different series of solutions are alternatively 
very close to being homogeneous. This can be seen in the two upper panels of the figure: For most of the represented 
values of L, there exists a solution whose energy is close to zero. The spatial profiles of such solutions reveal they are 
nearly homogeneous in /. 

Let us first consider the case of a positive detuning; see the left panels of the figure. For L « 0, there are two 
stationary solutions. The one with highest energy is analogous to the so-called first sol-sol solution in the fine-tuned 
case, in that it contains fractions of solitons attached at a; = —L and x = +L. The one with lowest energy is close to 
being homogeneous. When increasing L the “sol-sol” solution moves towards the other one, as the fraction of soliton 
at x = ±L decreases. Increasing L, we see a transition similar to an avoided crossing in quantum mechanics, and the 
series of solutions corresponding to the “sol-sol” becomes nearly homogeneous while the other solution shows fractions 
of shadow solitons at x = ±L, therefore becoming analogous to the “sh-sh” solution. Two new series of solutions with 
the same energy and related by a parity transformation appear at a critical value of L close to 1 for the parameters 









11 



FIG. 8. Energies of the first nonlinear stationary solutions of the GPE as functions of the half-distance L between the two 
discontinuities of the potential, for four different detuned sets of parameters. The parameters J = /(,,ext = 1, /p,ext = 0.7, 

and /i,,int = 1.5 take the same values for these four plots. The two plots on the left correspond to a positive detuning, with 
/p.int = 0.99 (top) and 0.9 (bottom). The two plots on the right correspond to a negative detuning, with /p,int = 1.01 (top) 
and 1.05 (bottom). 


of the figure. These series correspond to “sh-sol” and “sol-sh” solutions, in that both of them have part of a soliton 
at one end of the internal region and part of a shadow soliton at the other end. With increasing L again, the nearly 
homogeneous solution goes to the second “sh-sh” solution. We find the same pattern repeats itself periodically: One 
branch of solutions appears corresponding to a “sol-sol” one, which becomes more homogeneous, generates one series 
of “sh-sol” and “sol-sh” solutions, and turns continuously to a “sh-sh” solution. The set of stationary solutions is 
thus very similar to the fine-tuned case, except that no series of solution goes continuously from one “sh-sh” to the 
next “sol-sol.” Instead, an avoided crossing separates the two series. 

The case of a negative detuning is very similar, except that the avoided crossing is replaced by an anticrossing: The 
two initial solutions merge at a critical value of L and two new solutions appear at a second, larger critical value. 
Interestingly, even when considering higher-energy solutions not represented in the figure, no stationary solution exists 
between the first two critical values of L. In that case, numerical simulations using the same code as those presented 
in the next section always show an emission of infinite soliton trains. We thus recover a situation similar to the one 
found in Ref. [21], in that soliton trains arise from the absence of stationary solution. 

In brief, even though a small detuning introduces some modifications of the linear series of stationary solutions, such 
as introducing some avoided crossing or antiavoided crossing, it does not significantly affect the physical properties of 
the set of stationary solutions. 


IV. TIME EVOLUTION OF BLACK HOLE LASERS 

In this section, we study the time evolution of the density in black hole laser configurations. We first analyze 
individual histories associated with specific initial conditions. In order to relate these histories to the ensemble 
averaged density observed in Ref. [1], we study the evolution of the mean value over several solutions with different 
initial conditions. 

For definiteness, we present numerical simulations with tuned parameters. We explicitly checked that a small 
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FIG. 9. (Color online) Plot of / as a function of x for different times: t = 0 (top, left), 25 (top, right), 50 (bottom, left), and 
75 (bottom, right). The parameters are J = /p,mt = /t.ext = 1, /p,ext = 0.7, fb;mt = 1-5, and L = Lo + ^ ~ 0.68. The 

initial conditions are such that f{t = 0) — 1 > 0 for — L < x < L. This is similar to introducing a positive detuning. Notice 
that the range of <5/ = / — 1 in the first two plots is [—0.001,0.007], whereas it is [—0.02,0.25] for the last two. 

detuning does not change the main results. All numerical simulations presented below have been done on a large 
torus of length 4807r/u, where v is the velocity of the flow in the homogeneous solution. The integration was done 
on a uniform grid with 8196 space points and a time step of 5 x 10“^. We checked that dividing the space step by 
2 and the time step by 4 did not change the numerical solutions in a noticeable way. The initial conditions consist 
in a superposition of two waves of constant amplitude Sf/f = 10“^, and with a wave vector k = which 

is the dispersive zero-frequency root in the initially homogeneous supersonic region. These give a significant initial 
amplitude of the most unstable laser mode, and a relatively small amplitude to the real-frequency modes. We adopted 
these initial conditions for practical convenience, and we checked that similar results are obtained when using different 
initial conditions. 


A. Nonlinear effects on individnal confignrations 

To start, we ran simulations with L < Lq. As expected, the amplitude of the perturbation remained of order 10“^. 
As expected as well, non-trivial aspects appear for L > Lq. Figure 9 shows the case where Lq < L = Lq + Xo/8 < Li, 
and when f{x « 0, t = 0) — 1 > 0. This sign acts as a positive detuning in the sense that it sends the solution towards 
the stable “sh-sh solution” with n = 1. At early times® (though larger than l/Thom., where Thom, is the imaginary 
part of the complex frequency of the ABM on top of the homogeneous solution), the evolution is dominated by the 
laser mode which dictates both the shape of 6f{x,t) = f{x,t) — 1, and its exponential growth. We verified that the 
growth rate is equal to Thom. — 0.27, computed by solving Eq. (Bll). At later times, of order t ~ 30, one enters a 
nonlinear regime. The growth rate decreases and the shape of Sf is slightly modified as the solution is approaching 


® We use “early times” for times large enough for the unstable mode with the largest growth rate to dominate, but small enough for 
nonlinear effects to be negligible and “late times” for times large enough for nonlinear effects to be important. In contrast, in Ref. [1] 
“late times” refers to times where the most unstable mode dominates. 

























13 


f f 






FIG. 10. (Color online) Plot of / as a function of x for different times: t = Q (top, left), 25 (top, right) 50 (bottom, left), and 
75 (bottom, right). The parameters are J = a/8/3, /p,int = /s.ext = 1, /p,ext = 0.7, fb,mt = 1.5, and L = Lq + ^ « 0.68. 
The initial conditions are such that f{x, t = 0) — 1 has the opposite value as that of the former plot. Here we thus have the 
equivalent of an initial negative detuning in the internal region. One verifies that the final profile, obtained after having emitted 
the soliton, is identical to that of the former plot. 


the “sh-sh” solution. This process is smooth, in the sense that no large-amplitude perturbation is emitted away from 
the horizons, and Sf(0,t) is a monotonically growing function of time. At late time the flow is stationary, and the 
spatial profile is exactly given by that of the “sh-sh” solution; i.e., it is fixed by L, and not by initial conditions. 

Figure 10 shows the evolution for the same value oi L = Lq + Ao/8 when the initial value of 6f{x,t = 0) has the 
opposite sign. In this case one has the equivalent of a negative detuning. At early times, \6f\ grows exponentially 
with the same rate Fhom., exactly as predicted by the Bogoliubov-de Gennes equation. However, instead of making 
it saturate, nonlinear effects now turn the hollow of the “sol-sol” solution with n = 1 into a soliton which is emitted 
towards x oo. The remaining value of Sf(x « 0) is now positive and saturates on the same solution as above: 
the “sh-sh” solution with n = 1. Therefore, at very late times, the two solutions obtained by flipping the sign of the 
initial value of Sf both asymptote to the “sh-sh” solution, and this despite their different behaviors at intermediate 
times. Moreover, we verified that the convergence is exponential with a decay rate given by the imaginary part of the 
frequency of the quasi-normal mode (QNM) [6] on this solution. We performed simulations with other different initial 
conditions, and always found the same end state, which is thus an attractor. It should also be noticed that the Z 2 
symmetry, which is present at early times, is thus completely broken at late times. This has important consequences 
on observables which are odd in Sf, as shall be shown in the next subsection. 

We found similar results when choosing Li < L < L 2 . In this case, we also found that the end state corresponds to 
the ground state, the “sh-sh” solution with n = 1. At early times Sf grows exponentially, with a sign which depends 
on the initial conditions. The difference with respect to the previous case is that the frequency of the laser mode now 
has a non-vanishing real part, so that Sf periodically changes sign in the linear regime. As a result, the sign of Sf for 
t = 0 in the internal region is no longer directly related to the emission of a soliton. Importantly, we here observe the 
first manifestation of a general tendency. When increasing the interhorizon distance 2L, the set of lasing modes gets 
larger. Consequently, the behavior of nonlinear solutions becomes more intricate, and less straightforwardly related 
to the initial conditions. It is possible that this complexity will lead to a chaotic, i.e., unpredictable, behavior when 
there are several lasing modes. It would be extremely interesting to validate this conjecture. 

To verify that complexity increases with the number of laser modes, we ran simulations with L 2 n < L < T 2 (n+i) 
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FIG. 11. (Color online) Plot of / as a function of x for different times: t = 0 (top, solid), 50 (top, dashed), 100 (bottom, 
solid), 200 (bottom, dashed), 300 (bottom, dotted), and 400 (bottom, dot-dashed). To ease the reading, the last three curves 
on the bottom plot have been displaced upwards by 1 (dashed), 2 (dotted), and 3 (dot-dashed). The parameters are J = ^/S/S, 
/p.int = /i),ext = 1, /p,ext = 0.7, /i.int = 1-5, and L = Lq + |Ao ~ 3.0. At early times, for t = 50, in the internal region \x\ < 3 
one observes the laser mode with the highest growth rate. At intermediate times, one observes some solitons which are not 
equally spaced due to a transient behavior. At late times, they tend to be equally spaced. 


for n = 1, n = 2, and n = 6. In the three cases, at early times, we found that the laser mode whose frequency has 
the largest imaginary part dominates. As a result, f{x,t) goes very close to the mth “sh-sh” solution for times 
of a few l/Tm- Because this stationary solution is not stable, after a while, the time-dependent solution emits one or 
several solitons which escape to x —>■ oo and then approaches the Zth “sh-sh” solution, with I < m. In our simulations, 
we observed that the solution quickly evolves, and saturates close to the “sh-sh” solution with n = 2. (At present we 
do not clearly understand this observation. We suspect it is due to the fact that mode mixing across the horizons is 
large, as the transition is sharp since we are using discontinuous parameters. We thus conjecture that the solution 
will evolve more slowly when using smooth profiles for V and g.) After having approached the “sh-sh” solution with 
n = 2, we observe in Fig. 11 the emission of solitons in an apparently periodic way, as was already observed in Fig. 5. 
At present we have not been able to identify any criterion able to distinguish the solutions that shall emit soliton 
trains, from those which shall not.® Let us here note that similar soliton trains have been observed in Ref. [21]. 

In any case, these complex behaviors result from an interesting interplay between the linear instabilities governing 
early time behaviors and nonlinear effects at later times. While linear instabilities trigger the cascading between the 
various “sh-sh” solutions with decreasing values of n, going from one solution to the next one is always accompanied 
by the emission of solitons. Depending of the solution, either a finite number of solitons is emitted, or an infinite 
number of solitons, as for white hole flows, so that a stationary solution is apparently never reached. A categorization 
of the set of possible behaviors, and their respective domains in parameter space, is perhaps possible but beyond the 
scope of this work. 


® Carusotto, Finazzi, and de Nova [22] observed emission of solution trains in their numerical simulations of black hole lasers. Later, de 
Nova informed us that he had numerically found that (for a fixed initial uniform density /) there exists a L-dependent threshold value 
of fb,int ~ /p.ext above which infinite soliton trains are emitted. This is in agreement, and completes, our own findings. More work is 
necessary to identify the mechanisms which determine the threshold. 
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FIG. 12. (Color online) Time evolution of the ensemble average of 5p (with respect to the homogeneous solution) for the two 
simulations of Figs 9 and 10 (left), and for two simulations with the same parameters except that L = Lq + |Ao (right) so as 
to have a nondegenerate dynamical instability. The solid (blue) lines represent the mean (5p), whereas the dashed (orange) 
/ _ 2 \ 1/2 

lines represent the rms ( ((5p) ) . A bar means space average over the internal domain —L < x < L, and () means average 

over the two simulations with opposite perturbations at t = 0. The dotted (red) lines show exponentials with growth rates F 
and 2r, where F is the imaginary part of the frequency of the laser mode. One clearly sees that the mean grows with a rate 
which is twice that of the rms value. One also sees that the two quantities coincides at late times. 



FIG. 13. (Color online) Relative differences of the space averaged density perturbation, as functions of time. The left (respec¬ 
tively, right) panel corresponds to the left (respectively, right) panel of Fig. 12. The dotted (red) lines show exponentials with 
growth rates equal to the imaginary part of frequency of the QNM which lives on top of the ground state. At early times, the 
fiat plateau of height equal to 2 reveals that both perturbations Spi, i = 1,2 grow at the same rate, and stay equal and opposite 
to each other. It ends when nonlinearities become signihcant. At late time, the difference of Spi exponentially decreases, as 
both configurations approach the same ground state. Hence this difference is governed by the time dependence of the QNM on 
that ground state. 


B. Time evolution of the mean density, breaking the Z2 symmetry 

As mentioned in the introduction, Steinhauer observed that the averaged value of the density (taken over 80 
realizations) develops a clear spatial pattern with a rapidly growing amplitude; see Fig. 2 in Ref. [1]. The nodes of 
the profiles are fixed and compatible with those of the most unstable lasing mode, which dominates the growth of 
the 52 ; see Fig. 4 in Ref. [1]. In the following we show that a behavior similar to that of his Fig. 2 can be obtained 
from the breakdown of the Z 2 symmetry by non linear effects. A precise definition of this symmetry can be found in 
Appendix A 1 . 

To illustrate the roles of the Z 2 symmetry and its breaking, we show in Fig. 12 the space and ensemble average 
of ( 5/9 = p — 1 = /^ — 1, over the internal region —L < x < L. Our ensemble is very simple as it only contains two 
realizations with opposite initial perturbations, i.e., related by the Z 2 symmetry. We denote (•) the ensemble average 
over the two realizations, and ^ the space average over the region —L < x < L. In Fig. 12, the solid lines represent 

((5p). For comparison, we also represent the root mean square (rms) value \ \5p ) (dashed) whose behavior is closely 
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related to that of the 52 ; see Appendix A 2. These plots are obtained in flows possessing one degenerate instability 
(left) and one nondegenerate instability (right). 

At early times, we notice that (Sp) grows as This is due to the suppression of {6p) by the Z 2 symmetry 

to linear order: {Sp') is of order 2 in the amplitude of the perturbation, while the rms, which is not suppressed by 
the symmetry, remains linear in this amplitude. At late times instead, for both the left and the right panels, the two 
curves giving the mean and the rms values are indistinguishable. This is a direct consequence of the fact that the 

profiles of f{x), and thus those of p{x), of the two simulations become identical, so that (6p) = (Sp') = Sp , where in 
the last expression Sp is the common value of the density perturbation for the two simulations. In the nondegenerate 
case, because of the real part of the frequency of the lasing mode, the two density perturbations periodically vanish 
at linear order, hence the hollows in the plot of the rms. 

A complementary view of this symmetry breaking is shown in Fig. 13. In this figure we show the relative difference 
between the space averaged values of the density perturbations Spi as a function of time. At early times, its value 
remains close to 2, as the density perturbations remain opposite to each other while increasing exponentially. Its value 
deviates from 2 when nonlinear effects become important, and exponentially decreases at late times. We verihed that 
the late time decay rate is equal to the imaginary part of the frequency of the QNM defined on top of the “sh-sh” 
solution. 

In brief, in this subsection, we saw that the ensemble averaged value of the density fluctuation Sp exponential grows, 
even when its initial value is zero. The departure from zero is due to the breaking of the Z 2 symmetry by nonlinear 
effects. In this case, its growth rate is twice larger than that of the rms value of Sp, which is also that of Spi of each 
particular realization. The doubling of the growth rate of the mean value can thus serve as an unambiguous test to 
determine whether the non-vanishing value of Sp is due to nonlinear effects or to classical initial conditions. 


V. CONCLUSIONS 

The aim of this paper was to characterize the time-dependent and the nonlinear effects at play in the evolution of 
one-dimensional black hole laser flows. To this end, we first numerically studied the time evolution of transonic flows 
containing a single horizon. For simplicity, and to make contact with the analytically computed set of stationary 
flows, we worked in the steep-horizon limit. 

Even though the set of stationary solutions is the same for flows mimicking a black hole or a white hole, we saw 
that their time evolutions radically differ. Black hole flows evolve in such a way that the local values of the frequency 
w and current J are changed to reach a member of the linear series of stationary solutions which are asymptotically 
homogeneous on both sides. Even though we reached this conclusion in the steep-horizon limit, we believe that this 
result should also apply to potentials which have a smooth spatial profile. We checked this numerically on a few 
examples of smooth potentials. Consequently, we strongly conjecture that the stationary flow which is homogeneous 
(in that it contains neither undulation in the supersonic region nor soliton on the subsonic side) and has a value 
of the current compatible with the conservation laws acts as an attractor for transonic flows solutions of the GPE 
which mimic black holes. A first tentative proof was presented. If this could be demonstrated, it would mean that 
one-dimensional analog black hole flows also obey a no-hair theorem. 

Considering time-dependent white hole flows, we observed that they always emit an undulation in the supersonic 
region. Depending on the sign of the detuning (with respect to the asymptotically homogeneous solution on both 
sides) fixed by initial conditions, the undulation amplitude either saturates to a constant value when a shadow soliton 
is attached to the horizon, or widely varies, signaling that deep solitons are emitted in the subsonic region. This 
second scenario is found for a “negative” detuning, which sends the solution to the unstable soliton solution. 

When turning to black hole laser configurations, and when the initial density profile is smooth, at early times, the 
evolution is dominated by the most unstable laser mode. In this regime, time-dependent solutions obtained by flipping 
the sign of the initial density perturbations remain equal and opposite, until nonlinearities become significant and 
break the Z 2 symmetry. At late times, we observe a wide variety of behaviors. When the number of unstable modes is 
small, we found that the solution reaches the stationary ground state either smoothly or by emitting a finite number 
of solitons. When there are many unstable modes, we found that apparently infinite soliton trains are emitted. 

To make contact with the observations of Ref. [1], we then computed ensemble averaged observables. Because of 
nonlinearities, the mean density can acquire a non-vanishing value even when its initial value identically vanishes. 
This could provide an explanation for what was observed without relying on a bias in the initial conditions. Moreover, 
this scenario should be distinguishable from that based on initial conditions because in the former the growth rate of 
the mean is twice that of each particular realization. 

To go further, it would be interesting to precisely compare the results of Ref. [1] with our predictions to determine 
whether the breaking of the Z 2 symmetry comes from nonlinear effects, or from initial conditions. From a more 
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theoretical point of view, a study of soliton trains and their generation using analytical techniques would be very 
useful in determining the set of initial conditions leading to their formation. This would also tell us to what extent 
the evolution is predictable or chaotic in character, i.e., what can be said of the late-time evolution given initial 
conditions. Finally, it would be enlightening to see how our results generalize to systems with subluminal dispersion 
relations, such as surface water waves [9, 23]. As the key elements are unchanged, in particular the fact that an 
undulation cannot be emitted by transcritical flows mimicking black holes, we expect that a kind of no hair theorem 
will also apply to such flows. In addition, black hole laser configurations, which are now obtained by enclosing a 
subcritical region between two supercritical ones, should also behave in a manner similar to that analyzed in this 
work. 
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Appendix A: The Z2 symmetry and the behavior of the 32 correlation function 

1 . The Z2 symmetry 

As the Bogoliubov-de Gennes equation contains only linear and antilinear terms, its set of solutions is invariant 
under multiplication by —1. Moreover, this operation does not change the physical properties of the perturbation, 
such as its energy or momentum, as it amounts to a change of phase by tt. As a relative perturbation of the condensate 
wave function (p gives a relative density perturbation df/f = ^4>, we thus have a Z 2 symmetry Sf —)■ —Sf of Eq. (3), 
leaving the energy of the solution unchanged to linear order. Therefore, when working with a thermal state (or any 
other state which does not break this symmetry) the average of Sf, or of any observable which is odd in Sf, is and 
remains identically equal to zero. This applies both to the ensemble average value of the undulation amplitude emitted 
by white hole flows and to density fluctuations associated with the black hole laser instability. This Z 2 symmetry will 
be broken by nonlinear terms in the GPE, the first of which are quadratic in Sf. So, the ensemble average value of 
Sf will generally develop an expectation value of order <5/^. 


2 . Time evolution of the correlation functions 

We define the G 2 correlation function by 

G 2 {x,x';t) = {p{x,t)p{x,t)) - {p{x,t)) {p{x,t)), (Al) 

where {p{x, t)) is the average of the local density over a given set of realizations with different initial conditions. When 
divided by the product {p{x,t)) {p{x',t)) it gives the standard g 2 function. In our units, the mean density {p(x,t)) 
is always close to 1. Hence to a very good approximation, our G 2 {x,x'-,t) basically agrees with the standard (/ 2 - We 
chose to represent G 2 as the link with the time evolution of Sec. IV is more straightforward. 

When using ensembles made of only two realizations with opposite initial conditions on Sf, one clearly sees the 
effects, and the breakdown, of the Z 2 symmetry. We here use the same numerical simulations as those of Subsec. IV A, 
for Li < L < L 2 , i.e., in the case where there is only one nondegenerate laser mode. Figure 14 shows G 2 {x, x'; t) at two 
different times. At early times, we observe a pattern with two peaks and two hollows between the two horizons, which 
closely corresponds to that engendered by the dynamically unstable mode. (The left panel shows a slight asymmetry 
between the black- and white hole horizons, which seems to be due to a residual effect from the initial conditions.) 
By varying the time (not represented), we verified that the growth rate is equal to 2r, where F is the imaginary part 
of the frequency of the dynamically unstable mode. The right panel shows G 2 (x,x';t) at a time when saturation 
effects become important. At that time, in each of the realizations, the solution has moved close to the stable sh-sh 
stationary solution. The two main differences with respect to the previous plot are the profile in the internal region 
which is now typical of the sh-sh solution and the relatively important correlations between the regions x < —L and 
—L < X < L. These are due to the emission of phonons from the black hole horizon, i.e., to the Hawking effect [5]. 
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FIG. 14. (Color online) Correlation function G 2 {x,x';t) of Eq. (Al) for the parameters of Figures 9 and 10, except for 
L — Lo + |Ao ~ 1.3, at two different times: t = 20 (left) and 80 (right). 



FIG. 15. (Color online) Correlation function G 2 {x,x';t) of Eq. (Al) for the parameters of Figures 9 and 10, except for 
L = Lo + |Ao ~ 1.9, at t = 10 (left) and t = 200 (right). The sharp features observed for x,x' ~ A on the right panel are due 
to the emission of a soliton. 


When the two solutions saturate, they become nearly identical and the G 2 becomes accordingly very small, with a 
typical amplitude equal to the square of that of the real-frequency modes present in the initial conditions. In fact 
when t = 120 we found that the peak-to-peak value of G 2 is less than 10“^. 

Two relatively important differences arise when more unstable modes are present. First, the early-time behavior 
is more complicated as the configuration of the left panel of Fig. 14 is reached after a series of steps progressively 
reducing the number of nodes in the supersonic region, as explained in Subsection IV A. An example of such an 
additional step is shown in the left panel of Fig. 15. Figure 15 is obtained for L « 3.0, so that two dynamically 
unstable modes are present. On the left panel, one sees that the density-density correlation function shows a different 
pattern with five peaks and four hollows. Second, at late times one or the two solutions may generate a seemingly 
infinite soliton train. In that case G 2 {x,x';t) remains large at late times; see the right panel of Fig. 15. The peak for 
a; « x' « 5 is due to a soliton emitted in the internal region and expelled to infinity, whereas the peak at x « x' « 0 
is due to a soliton which is being produced in the internal region. 

In brief, we see that the analysis of the time-dependence of the density-density correlation function reveals the 
various steps of the dynamical evolution of black hole lasers, from the growth of the most unstable mode at early 
times to the various scenarios at late times. 


Appendix B: Stationary solutions and ABM 

In this appendix we explain how the stationary solutions and ABM can be found in black hole laser configurations 
in the steep horizon limit. 
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1 . Stationary solutions in the detuned case 


Following Ref. [6], we find the stationary solutions in each of the three regions Ii : x < —L, I 2 ■ —L < x < L, and 
I 3 : X > L, and match them at a: = ±L to find the solutions on K. We also assume that the flow is uniform and 
subsonic for x —>■ ± 00 . In each region Eq. (4) gives 

= + |r-^+C., (Bl) 

where Ci is an integration constant. Imposing that the solution be homogeneous and subsonic in the limits x —> ±00 
gives Cl = C 3 = C„iax,ext- On the other hand, in the internal region C 2 can be varied continuously. The phase 
portrait of Eq. (4) is shown schematically in Fig. 16 for different values of the parameters gint, Mint, and Cnt- The 
red curve shows the trajectory in phase space (/, p = dxf) of the solution in the external regions, while the blue one 
shows its trajectory in the internal region. Point A corresponds to the homogeneous supersonic solution in the internal 
region and point C to the homogeneous subsonic solution in the external regions. So, the tuned case corresponds to 
A = C. Global solutions are found by following the red line in the direction of the arrows from point C to one of the 
intersection points with the blue line (R,Z1, if,E), then to another one or the same intersection point following the 
blue line, and then back to C following the red line. The first step corresponds to the region Ii, the second step to 
I 2 , and the third one to I 3 . For each solution, the required length of the internal region is given by 

2T= (B2) 

J P 

the integral being evaluated over the path followed in the second step. The procedure is the same as that used in 
Ref. [6], to which we refer for more details. The main difference in the presence of detuning is that there is no 
homogeneous solution. The set of solutions thus qualitatively changes for values of Cint at which the number of times 
the blue curve crosses the red one changes. To express these critical values, it is convenient to first define 


/*s,ext — 


24%xt 


f2 " I i2 ■ 
•/6,ext ' Jp,ext 


(B3) 


/s.ext is the value of / at the bottom of the stationary soliton in the external regions. The first critical value of Cint 
is the one for which the blue line is tangent to the red one at point E = B. It is given by 


Cn 


= Ce: 


(Mint Mext) 

2 (Pint Mext) 


(B4) 


C'int.m is the minimum value of Cint for which the matching conditions, i.e., continuity of / and dxf at x = ±L, 
can be satisfied. In the fine-tuned case, it is equal to Cnt,min- The second critical value of Cnt is the one for which 
E = F = C (for a positive detuning) 01 B = D = C (for a negative detuning). It is given by 
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(/^int + fp.int + flintflint) ^ (Z^int + /p.int) 
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Another critical value is the one for which the blue and red lines are tangent at B = D, and is given by 

Cnt,s = Mext/s,ext “ ^!^/s,ext + 


2/Z 


s,ext 


(B5) 


(B6) 


Filia,lly, the last critical value of C^lnt is C^int,max- Fi§- 16 shows the four cases C*int,m C^int C^int ,0 5 ^int,S 5 C^int,max 
(two upper panels), (Fint,m ^ C^int,0 C^int ^ C^int,s; C^int,max (bottOm left panel), and Cint,m C^int,0 7 ^int ^ 

C'int,max (bottoui right panel). 

The energy of a stationary solution is defined through 


^ = J dx + Eq, 


(B7) 


where Eg is a constant. Using the GPE, and adjusting the constant Eq so that the energy of a uniform solution with 
density ext vanishes, this may be written as 


E = -j ^-^{f{xf-fUt)dx. (B8) 

We found no important difference when using the other energy functional of [6], defined by a Legendre transform 
of Eq. (B7) to impose the value of the current J. 
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FIG. 16. (Color online) Schematic drawings of the phase portraits p = dxf vs / of Eq. (4), restricted to C = Cext.max in the 
external regions (red line) and one value of C between Gint.min and Cint.max in the internal regions. Arrows give the direction 
in phase space for increasing values of x. The 4 panels show the trajectories in phase space for four different sets of parameters. 
The detuning manifests itself in the separation between the points A (located at / = fp,int,P ~ 0) and C (/ = /f),ext,P = 0). 


2 . Finding the complex-frequency modes 

We remind the dispersion relation of perturbations in a region with uniform density p and flow velocity v is 

(A - vkf = gpk'^ + ^, (B9) 

where A is the angular frequency and k is the wave vector. When A € C — K, the four solutions in k of Eq. (B9) 
are complex. It is easily shown that two of them have a positive imaginary part while the other two have a negative 
imaginary part. Let us denote as ki, i G {1,2, 3, 4} these four roots, with < 3^2 < Qfca < $ 5^4 and as (j)i the 
solution of the Bogoliubov-de Gennes equation in a homogeneous background with wave vector ki. In our black hole 
laser configuration, for x —>■ —00, (f>i and (j )2 decay exponentially while ^3 and (/)4 grow exponentially. On the other 
hand, for x -G 00, 4>^ and tpi decay exponentially while (j)i and (p 2 grow exponentially. Let be the solution of 
the Bogoliubov-de Gennes equation which goes to cj)i for x —>■ ±00. We define the coefficients Aij, {i,j) G (1,2,3,4}^ 
through 

4 

4>.,_ A,,. (BIO) 

1=1 

An acceptable solution must be a linear superposition of $1 _ and ^ 2 - (to be asymptotically bounded at x —)■ —00), 
and a linear superposition of and 4 > 4 __|_ (to be asymptotically bounded at x —>■ -boo). Such a solution exists if 

















and only if the determinant 
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^ 1,1 

^24 


^ 1,2 

^ 2,2 


= 0 . 


(Bll) 


Our strategy is thus to look for zeros of the determinant in the complex plane. To this end, we found it convenient to 
compute its phase along closed lines. Indeed, this determinant is holomorphic in the complex plane, except for branch 
cuts which are easily identified and terminate at values of A for which the dispersion relation has a double root. By 
choosing a contour which does not cross a branch cut, the phase shift is equal to 27rn, where n is the number of zeros 
of the determinant inside the contour, which can then be refined to locate them more precisely. In practice, we chose 
a rectangle with vertices at iwmax and iwmax + * Tc, where 

a;inax = |u| + + 8cM ^ / 2 , q 2 ) ’ 

V3 |u| + + 8c^ / 


evaluated at a; = 0, and Tc is the imaginary part of the same quantity evaluated at a; —>■ ± 00 . The frequencies of 
all the asymptotically bounded modes we found are well inside this contour, and we saw no evidence of complex 
frequencies outside while looking at the evolution of the phase along the contour (as a nearby zero would give a rapid 
variation of the phase) or by extending the contour. When dynamical instabilities were found, the contour was then 
refined to locate them with an accuracy of 10% for both the real part and the imaginary part. 

The method used in Ref. [6] worked only for perturbations of the homogeneous solution. It was found that its degree 
of stability changed each time a new “connected” series of solutions with negative energy appears. The corresponding 
critical values of L are 
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Even values of m correspond to the appearance of a degenerate dynamical instability, while odd values correspond to 
a degenerate instability turning into a nondegenerate one. Figure 6 shows the degree of stability of the “connected” 
stationary solutions. We first remark that apart from the homogeneous solution for L < Lq, the only connected 
dynamically stable solution is the first “sh-sh” solution of Ref. [6], i.e., the one with lowest energy. We checked 
numerically that the “non-connected” solutions are all dynamically unstable. The first “sh-sh” solution is thus the 
only dynamically stable stationary solution. This partially proves the conjecture that was formulated in Ref. [6], 
namely that if the system evolves towards a “connected” stationary solution at late times, then the final state is the 
solution with lowest energy, i.e., the homogeneous solution for L < Lq and the first “sh-sh” solution for L > Lq. The 
second part of the conjecture, namely that the system generally becomes stationary at late times, is investigated in 
Sec. IV. The first “sol-sol”, “sol-sh”, and “sh-sol” solutions all have a degenerate dynamical instability. We found 
that, in general, a nondegenerate dynamical instability appears when going from one “sh-sh” (respectively, “sol-sol”, 
“sh-sol”, or “sol-sh”) solution to the next one. This could be expected from the results of Ref. [6], where it was 
shown that the homogeneous solution gains one nondegenerate dynamical instability each time L is increased by 
Ao/2. Our present numerical calculations confirm that series of solutions which can be continuously deformed into the 
homogeneous one inherit these additional instabilities. The only exception we found is the second series of “sol-sol” 
solutions, as for some values of L two solutions of this series coexist. Then, as shown in the insert of Fig. 6, the one 
with lowest energy has a degenerate and a nondegenerate dynamical instabilities, while the one with highest energy 
only has a nondegenerate instability. The same pattern repeats itself for the next series of “sol-sol” solutions, with 
the addition of one nondegenerate dynamical instability when going from one series to the next one. 

A small detuning has little effect on the set of stationary solutions, except for L « Lm, m gN. Stationary solutions 
can thus be identified with the ones studied above. We found that the above results on linear stability continue to 
hold; see Fig. 8. We also conjecture that they remain true for smooth variations of g and /r. Although such setups 
could in principle be examined using the method described above, we leave this to a later work. 

To end this section, let us compare the growth rate of the unstable modes. In Ref. [6] it was shown that, unless 
L is very close to one of its critical values, the most unstable ABM is the one which appears last, its growing rate 














22 


being larger than that of the second most unstable mode by a factor of order 10. We found a similar result for the 
non-homogeneous solutions. Moreover, when considering an inhomogeneous solution which appears for L = L 2 rm 
m S N, we find that the set of complex frequencies on this solution is close to that on the homogeneous solution for L 
slightly below Lm- In other words, the new series of solutions inherits the ABM that were present on the homogeneous 
solution for L < Lm- It is thus less unstable than the homogeneous solution, which has a new dynamical instability 
with a generally larger growth rate. This is different for series of solutions appearing at L = L 2 m+i, as when crossing 
this critical value of L no new unstable mode appears on the homogeneous solution. Instead, a degenerate instability 
is converted to a nondegenerate one, while on the “sh-sol” solution it remains degenerate. In that case we found the 
growth rates have the same order of magnitude, with the growth rate of the mode on the inhomogeneous solution 
being in general larger than that on the homogeneous solution; see Fig. 7. 


Appendix C: Some properties of dispersive shock waves 

In this appendix we determine the main properties of the dispersive shock waves emitted by a generic (detuned) 
initial black hole configuration. Our goal is to show how simple and general arguments explain the observations of 
Subsec. IIB. A more accurate study of dispersive shock waves and their generation can be found in Refs. [24, 25]. 
We notice on Fig. 3 that there is only one shock wave in the subsonic region. Other simulations we ran with different 
parameters also showed only one shock wave in this region. Knowing this, the velocity Vs of the shock wave and Vtun 
of the condensate behind it are easily computed using mass and momentum conservation. We find 


y |/|n (^1 + « ^^0±C-, (Cl) 

and 

Vfin= ^Vo+ (l- ^]v,, (C2) 

J fin \ fin / 

where vq (respectively /o) is the initial value of v (respectively /). The approximate equality in Eq. (Cl) is obtained 
for /fin ~ /o- To describe the shock wave emitted in the subsonic region, one must choose the solution with the minus 
sign. The new value of the frequency w^n may then be computed using the stationary CPE: 

2 (c^fin - K-) /fin - 2g_ /g\ - = 0, (C3) 

J fin 


where 


Tfin = ^'fin/fin• (C4) 

In the supersonic region, the two solutions of Eq. (Cl) are positive. So, we expect that two shock waves will be 
emitted. This is indeed what we saw in our simulations; see, for instance. Fig. 3. Using the same conserved quantities, 
we find that the density between the two shock waves, /fnt, is a solution of 


Vq ± 


ifS + fL) 


2/iit/o 


= Infill ± 


9+ifL-fL) (fin+ fL) 


2/iit/fin 


— ^int5 


(C5) 


where the two ± signs are independent and iffnt is the corresponding flow velocity. For a small detuning, we obtain 


/int « 2 

The velocities of the two shock waves are given by Eq. (Cl), with /o replaced by the value of / on the left of the 
shock wave, and /gn by the value of / on the right of the shock wave. 

Let us now consider the decay of the oscillating tail of the shock wave at late times. To this end, we first determine 
the evolution of a wave packet in a homogeneous background, and then impose reflexive boundary conditions at a; = 0. 
We write the condensate wave function as 


t) = e 


,i{—ujt-\-vx-\-59{x,t)) 


{fo + Sf{x,t)), 


(C7) 
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where is a solution of the GPE. To first order in Sf, its derivatives, and the derivatives of 59, the GPE 

gives 


Let us write 


{dt + vd^ f - gfodl + ]5f = 0. 


/ OO 

e’^^Ak{t)dk. 

-OO 


The time evolution of the functions Ak are given by 


{dt + ivkf + gf^k"^ + ^ ) Af,{t) = 0. 


The general solution is 




where are two complex numbers and 


jU4 


To be specific, let us start from the initial configuration 

5f{x,t = 0) = 


0 X < OV X > L 
Fq 0 < X < L 


(C8) 

(C9) 
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This is analogous to the detuned black- or white hole case in the limit L —)■ oo. Then, 
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At late times, we can do a stationary phase approximation. The terms in describe a wave emitted at the 

point X = L. We shall drop them as we are interested in the case where the limit L —>■ oo is taken before t —)■ oo. 

doj ^ 

In a stationary phase approximation, the only wave vectors which contribute are those for which ~ J- 
straightforward calculation gives the two possible solutions 


k = 



(C17) 


If the flow is subsonic, these wave vectors and the associated frequencies are imaginary in the limit x/t —)■ 0. The 
solution thus decreases exponentially. If the flow is supersonic, these solutions are real and we find in the limit 
x/t —>■ 0: 
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2tt kc (vkc — oJc) 
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Let US now explain how this calculation relates to the numerical observations. Its first prediction is that, at late 
times, the supersonic region should show a modulation with a wave vector kc- This is consistent with our observations; 
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see Fig. 3. Although we performed the calculation explicitly using very particular initial conditions, it is clear that 
this result is much more general. Indeed, as long as the initial conditions are such that A±k^ 7 ^ 0 and the stationary 
phase approximation holds at late times, the only modes which contribute for x/f —>■ 0 are those with a vanishing 
group velocity, i.e., with wave vectors zLkc and frequencies ±Wc- However, Eq. (CIS) predicts that the amplitude 
of the modulation does not depend on x, which is in contradiction with our observations. Moreover, according to 
Eq. (CIS) the amplitude decreases as Our numerical simulations in the black hole case show a stronger decay 

in The reason for these two differences is that Eq. (CIS) was derived assuming homogeneous potential and 

coupling constant. To mimic a black hole configuration, we must impose some boundary conditions at a: = 0 coming 
from the subsonic character of the flow for a; < 0. As Sf decays exponentially in this region, one can try to enforce 
reflective boundary conditions. We will see that such boundary conditions indeed give the correct result. Eor a; > 0 
with x/t^v, y/gfo, 5f is the sum of two waves: 

1 . the wave emitted directly from a: = 0 to a;, which has the form 

6 fi{x, t) = F{k) sin(fca; — + 7 r/ 4 ), (C19) 

where F(k) is the amplitude from the stationary phase approximation, evaluated at A: > 0 such that the group 
velocity is equal to x/t] 

2 . the reflected wave, which has the form 

Shix, t) = —F(k') sin(fca; — uj'^t + 7 r/ 4 ), (C20) 

where k' is the second solution of the dispersion relation close to for the same frequency lo. 

Here a few remarks are in order. First, the phases of the two waves are the same since the emission point coincides 
with the one where reflexive boundary conditions are imposed, so that the two waves propagate in the same medium 
and in the same direction at all times. If the reflexive boundary conditions were imposed at a point x = —e < 0, there 
would be a dephasing of (k — k') e. Second, the amplitudes of these waves are different since the amplitude given 
by the stationary phase approximation depends on k. Assuming again for a moment that the reflexive boundary 
conditions are imposed at a; = —e < 0, the incident wave at x = —e would have an amplitude F{k'), so the amplitude 
of the reflected wave is —F(k'). Assuming that F is differentiable at k = kc with a non-vanishing derivative, which 
can be easily checked by an explicit calculation, we thus get 

(5/(x, t) « (fc — k')F'{kc) sin(fcx — -I- 7 r/ 4 ). (C21) 

One can now extract the time and space dependence of the amplitude using that F' inherits the dependence in 
of F, while k — k' Ki 2{k — kc) « 2vg/^^, evaluated at fc = kc- Using that vq = x/t and that is finite and 
non-vanishing at k = kc thus gives 

(5/(x,t) cx sin(fcx - -I- 7 r/ 4 ). (C22) 

The phase 7 r /4 is not relevant as it depends on the precise initial configuration and boundary conditions. However, 
this argument tells us that, because k — kc scales as x/t, the sum of the incident and reflected waves gives a factor 
oc x/t to the amplitude of df with respect to the case of a homogeneous background on K. The result is, at late 
times and close to the black hole horizon, a sinusoidal perturbation of wave vector kc, angular frequency ojc and an 
amplitude scaling as This is fully consistent with the results of our numerical simulations. 
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